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NOMENCLATURE 


cp  heat  capacity 

£i  exponential  integral 

F(y)  defined  by  eq  4 

H  dimensionless  heat  flux  defined  by  eq  23 

k  thermal  conductivity 

f  1//3,  characteristic  length 

q  heat  flux 

R  defined  by  eq  33 

/  time 

T  temperature 

V  a  variable  defined  by  eq  10 

x  coordinate  in  direction  of  heat  flow 

y  defined  by  eq  2 

g  density 

a  thermal  diffusivity 

d  dimensionless  temperature  defined  as  7VTr 

V  defined  as  y/( 

t  dimensionless  time  defined  as  otot/H2 

t'  dimensionless  time  defined  by  eq  32 

b  thermal  boundary  thickness 

f  dimensionless  thermal  boundary  thickness  defined  as  6/f 

d  coefficient  in  the  expression  k(x)  =  fr0(l  +  fSx) 

7  dimensionless  distance  defined  as  x/(  or  x0;  or  Euler’s  constant 

Subscripts 

0,r  reference  state 

s  surface 


Approximate  Solutions  of  Heat  Conduction  in  a 
Medium  with  Variable  Properties 

YIN-CHAO  YEN 


INTRODUCTION 

A  number  of  books  deal  with  heat  conduction  in  various  geometries,  in  one  or  multiple  di¬ 
mensions,  and  with  constant  or  variable  properties  that  are  either  a  function  of  temperature 
or  position.  One  classic  work  that  deals  exclusively  with  conduction  is  Carslaw  and  Jaeger 
(1959),  in  which  all  the  problems  are  solved  analytically.  But  even  in  these  simple  problems, 
the  solutions  are  often  expressed  in  terms  of  integrals  and  series,  so  that  considerable  compu¬ 
tational  effort  is  needed  to  obtain  the  numerical  results. 

The  work  reported  here  represents  an  effort  in  obtaining  a  simple  but  reasonably  accurate 
method  for  the  prediction  of  temperature  distribution  in  a  medium  with  variable  properties. 
The  approximate  heat  balance  integral  method  (HBIM)  introduced  by  Goodman  (1958, 
1964),  which  has  shown  its  applicability  to  a  number  of  heat  conduction  problems,  is  extend¬ 
ed  to  the  present  case  of  a  medium  with  variable  properties.  The  concept  of  the  integral 
method  is  based  on  the  same  principles  involved  in  the  momentum  integral  method  of  the 
boundary  layer  in  fluid  mechanics  initially  proposed  by  Pohlhausen  (1921)  and  von  Karman 
(1921).  The  essence  of  this  method  is  that  at  any  time  t  the  thermal  boundary  layer  will  be 
limited  to  a  boundary  layer  thickness  6(1).  Beyond  this  depth,  the  temperature  will  remain  at 
its  initial  temperature  and  there  will  be  no  heat  transfer  beyond  this  point.  This  solution 
method  is  analogous  to  the  momentum  integral  method  in  that  the  basic  equations  are 
satisfied  on  the  average  over  the  volume  of  thickness  6(0  rather  than  at  each  point. 


MATHEMATICAL  ANALYSIS 
Constant  surface  temperature 

Consider  the  case  of  a  semi-infinite  solid  medium  initially  at  uniform  temperature  (let  it  be 
zero  for  convenience  in  the  mathematical  treatment)  with  its  physical  properties  given  as  a 
function  of  position.  For  1  >  0,  the  surface  is  maintained  at  a  temperature  of  Ts,  and  the  uni¬ 
directional  heat  conduction  equation  is  given  as 


d_ 

dx 


k(x) 


dT 

dx 


=  cpMeW  Yt 


(l) 


with  initial  and  boundary  conditions  as 


T  = 

0, 

for 

x  >  0 

and  t  s  0 

T  = 

Tv 

for 

o 

II 

* 

and  t  »  0 

•  1 


where  *  =  coordinate  in  the  direction  of  heat  flow 
k  =  thermal  conductivity 
cp  =  heat  capacity 
q  =  density 
t  =  time. 


A  new  dependent  variable  v  is  defined  as 


where  .v '  is  the  dummy  variable  and  kc  is  the  thermal  conductivity  at  some  arbitrary  refer¬ 
ence  point  (it  is  chosen  in  Mich  a  manner  that  the  subsequent  mathematical  treatment  will  be 
simplified).  By  substituting  eq  2  into  eq  1  and  simplifying,  eq  1  becomes 


<VT  _  F(y)  i)  T 
(iV:  nt0  tit 


where 


F(v)  = 


k(x)ev(x)  e(-v) 

CpO  Qo 


By  rewriting  eq  3  in  dimensionless  form,  this  becomes 


<VH  ;in 

H  ’  f""  T, 


Tr  and  f  are  the  reference  temperature  and  the  characteristic  length  of  the  medium  under 
study,  respectively.  A  new  dependent  variable  V  is  introduced  and  defined  by 


=  }F(y)0(y)  dy. 


By  using  eq  10  to  evaluate  the  derivatives  of  0  with  r;  and  r,  eq  6  can  be  expressed  in  terms  of 

V  by 
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d!V 
dr  drj 


(11) 


AM 

|  Hv) 


3j_K 

dr)1 


F'jr,)  dV 
F(r,) 


£|}- 


where  F'  =  dF/dr,.  The  corresponding  initial  and  boundary  conditions  are 


dv 

dy 


II 

© 

>3 

II 

o 

o' 

ll 

for 

r,  <  0, 

and  t  <  0 

^(0)^, 

for 

V  =  0 

and  t  >0 

where  3S  =  Ts/Tr  Integrating  eq  11  from  0  to  t  (which  is  defined  as  6/(  where  6  is  thermal 
boundary  thickness),  it  becomes 


d2V  | 

F'(e )  3T  | 

1  32T 

1  ^'(0) 

dV  ,  | 

dn2  ' 

,  F(e)  dr,  j 

,J  F(0)  dr, 2 

1  0  F(  0) 

3,  I  cj 

=  ^[K<«)-  k'(O)]. 


(12) 


Since  the  ..  ffect  of  surface  temperature  at  any  given  time  is  confined  to  r,  <  r,  it  is  expected 
that  the  following  conditions  will  be  satisfied  to  smooth  the  temperature  profile: 


II 

© 

•3 

II 

6 

(13a) 

33 

3rj  =  °’  "  = 

e 

(13b) 

Oj  I  Qj 

II 

O 

•3 

II 

f . 

(13c) 

The  corresponding  requirement  on  V  can  be  found  as 

V  =  0, 

3  =  0 

(14a) 

3  V 

T,  '  me- 

r)  =  0 

(14b) 

£  -  - 

T  =  f 

(14c) 

r)  =  f- 

(14d) 

f  can  he  represented 

as  a  cubic  polynomial  of  >/: 

1  a  •  hr,  »  Or 

*  d^' . 

(15) 

With  the  application  of  eq  14a  d,  the  coefficients  </,  h,  c  and  d  were  found  as 


a 


0 


( 1 6a I 


3 


■  « ■» 
-  ~  ■  t  „  *  ■  "  - ■ 


t.*!v 


Substituting  eq  16a-d  into  eq  15,  it  follows  that 


1  /■'«))«>, 


With  eq  17  and  14d,  eq  12  can  be  integrated  into 


F( 0)  ,  F(0)  1  ,  F( 0)  f 

riO)  (  m  in  - 2 - 


If  we  take  the  reference  properties  of  k,  cpand  q  to  be  those  at  the  surface,  then  HO)  =  1, 
and  eq  18  can  be  further  simplified  to 


2  ,  In  2  V— '  =  3f '(0)r. 


Ihe  computed  values  of  <  as  a  function  of  t,  with  F'( 0)  ranging  from  0.01  to  2.0,  are  shown 
in  ruble  1  and  in  Figure  1  for  selected  values  of  F'(0).  From  eq  10  we  have 

'2  F{'))H(v).  (2i 

Subsequently  from  eq  17  and  20,  the  temperature  profile  is 


»  =  +  i 

Hv)  \  f  ( 


(  (instant  surface  heat  flux 

I  he  same  analysis  as  in  the  case  of  constant  surface  temperature  can  be  applied  to  the  case 
of  constant  surface  heat  flux.  The  corresponding  boundary  conditions  on  1  become 


iM  ,  ill 

I  (0) 

m; 


//,  n  -  o 


where  H  is  the  dimensionless  heat  flux  defined  as 


in  which  q  is  the  heat  flux.  The  temperature  distribution  is  given  as 


/  <e)(2  *  /  (OH 


♦  JkrV.V.V'-  \  V 


Table  1.  Computed  values  of  dimensionless  boundary  layer  thickness  c 
as  functions  of  dimensionless  time  r  and  F'( 0)  for  a  constant  surface 
temperature  and  /TO)  =  1.0. 


_ FT 0) _ 

0.01  0.05  0.10  0.20  0.05  1.0 _ 2f) 


0 

0 

0.01 

0.348 

0.05 

0.775 

0.10 

1.102 

0.50 

2.460 

1.00 

3.479 

2.00 

4.940 

4.00 

6.760 

6.00 

8.605 

8.00 

9.960 

10.00 

11.150 

14.00 

13.250 

20.00 

15.900 

0 

0 

0.348 

0.348 

0.779 

0.785 

1.105 

1.115 

2.500 

2.551 

3.565 

3.667 

5.095 

5.310 

7.330 

7.752 

9.095 

9.726 

10.620 

11.460 

11.980 

13.040 

14.420 

15.901 

17.550 

19.730 

0 

0 

0.351 

0.357 

0.795 

0.825 

1.136 

1.198 

2.655 

2.973 

3.875 

4.529 

5.731 

7.075 

8.613 

11.390 

11.038 

15.295 

13.228 

18.995 

15.273 

22.585 

19.073 

29.510 

24.350 

29.560 

0 

0 

0.367 

0.388 

0.878 

0.987 

1.304 

1.528 

3.537 

4.750 

5.695 

8.222 

9.500 

14.758 

16.445 

27.348 

23.057 

39.710 

29.515 

51.975 

35.885 

64.188 

46.375 

88.500 

67.090 

124.850 

in  which  e  is  given  by  the  transen- 
dental  expression 

4  [1  _  1  1 

e  F 1  (0)1 2  2  -i-  ^ '  (0)e  | 

=  3F'(0)r.  (25) 

The  computed  values  of  t  as  function 
of  r  are  shown  in  Table  2.  Similarly 
a  plot  of  r  vs  e  is  shown  in  Figure  2 
for  some  selected  values  of  F'( 0). 
Appendix  A  shows  the  derivation  of 
eq  25. 

Comparison  with  exact  solution 

For  the  case  of  constant  surface 
temperature  ( T s  *  0)  that  was  ini¬ 
tially  at  zero  temperature,  the  exact 
solution  for  the  temperature  distri¬ 
bution  in  a  semi-infinite  solid,  with 
k(x)  =  Ar0(l  +0x)  but  with  q  and  cp 
being  constant,  has  been  reported 
by  Carslaw  and  Jaeger  (1959).  The 
numerical  solution  reported  by  Jae¬ 
ger  ( 1 956)  is  compared  with  those  ob¬ 
tained  by  this  approximation  method 


Figure  1.  t  vs  t  for  constant  surface  temperature, 
F  (0)  =  /. 


for  a  linear  dependency  of  fc(x)  on  x.  The  independent 


variable  y  becomes 


5  dx 

y  ~  Jd  +0x) 


j  ln(l  +0x). 


(26) 


5 


Table  2.  Computed  values  of  dimensionless  boundary  layer  thickness  t 
as  functions  of  dimensionless  time  r  and  F '  (0)  for  a  constant  surface 
heat  flux  and  F( 0)  =  1.0. 


F  (0) 


T 

0.01 

0.05 

0.10 

0.20 

0.05 

1.0 

2.0 

0 

0 

0 

0 

0 

0 

0 

0 

0.01 

0.246 

0.246 

0.246 

0.248 

0.253 

0.261 

0.277 

0.05 

0.549 

0.552 

0.555 

0.563 

0.587 

0.628 

0.718 

0.10 

0.775 

0.782 

0.790 

0.806 

0.853 

0.939 

1.131 

0.50 

1.740 

1.770 

1.809 

1.890 

2.148 

2.638 

3.792 

1.00 

2.465 

2.526 

2.607 

2.768 

3.313 

4.376 

6.875 

2.00 

3.495 

3.618 

3.778 

4.116 

5.275 

7.582 

12.930 

4.00 

4  959 

5.209 

5.554 

6.245 

8.745 

13.745 

24.965 

6.00 

6.092 

6.468 

6.968 

8.066 

12.005 

19.815 

36.980 

8.00 

7.198 

7.555 

8.233 

9.732 

15.170 

25.860 

48.980 

10.00 

7.898 

8.534 

9.390 

11.318 

18.285 

31.890 

60.990 

14.00 

9.379 

10.276 

11.505 

14.283 

24.430 

41.915 

84.980 

20.00 

11.259 

12.560 

14.358 

18.485 

33.581 

61.960 

120.990 

If  the  characteristic  length  is 
then 

i,=^-ln(l  +  >)  (27) 

or 

1  +  7  =  C  (28) 

where  ■>  =  x$  and  the  function  F( tj) 
becomes 

F(v)  =  (I  +  M  =  (!+*>)  =  C.  (29) 

For  the  case  of  k(x)  =  Ml  +  &x),  it 
then  follows  that  F(0)  =  \,F'(Q)  = 
1  and  F'(ri)  =  e’,  and  eq  21  reduces 
to 


for  the  case  of  constant  surface  tem¬ 
perature.  Similarly  for  constant  sur¬ 
face  heat  flux,  eq  24  becomes 


The  exact  solution,  in  this  case,  is  obtained  by  extending  the  solution  given  by  Jaeger  (1956) 
with  the  following  substitutions.  The  temperature  at  x  and  l  is  given  by 


in  which  R  =  r/a  where  a  is  the  radius  of  the  cylinder.  R  can  be  rewritten  as 
R  =  0+7)*. 

But  from  eq  27 
t) l  =  ln(l  +  7). 


(33) 


(34) 


Therefore  R  =  (e’O*. 

Since  q  and  cp  are  constant  and  can  be  considered  as  the  values  taken  at  the  surface,  the  r' 
defined  in  eq  32  can  be  expressed  in  terms  of  t  as 


t 

4(1 //j)2 


4  ' 


(35) 


The  tabulated  values  of  Jaeger  (1956)  for  t'  =0.1,  1 .0  and  5.0  (corresponding  to  the  t  val¬ 
ues  in  Jaeger’s  table;  corresponding  values  of  r  in  our  approximation  method  are  r  =  0.4, 
4.0  and  20.0  respectively)  are  shown  in  Figure  3  tc  ricilitate  the  computation  of  0/6s  for  var¬ 
ied  values  of  n  with  eq  34  to  get  the  analytical  results.  Figure  4  shows  comparisons  of  dimen¬ 
sionless  temperature  distribution  of  9/6s  vs  ??.  Only  the  case  of  constant  surface  temperature 
b  is  shown  for  r'  =  0.1,  1.0  and  5.0  (r  =  0.4,  4.0  and  20  respectively).  For  each  specific  t 
value,  a  corresponding  f  for  the  specific  value  of  F'(0)  =  1  is  obtained  from  Table  1  and 
subsequently  substituted  into  eq  30  to  evaluate  the  ratio  of  Q/Q%  as  a  function  of  >j.  In  Figure 


Figure  4.  Comparison  of  q  vs  6/6, 
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c.  T  =  20.0  (t'  =  5.0). 


Figure  4  (cont’d.). 


4,  it  can  be  seen  that  for  small  values  of  r  (i.e.  r'  =  0.1,  or  7  =  0.4),  the  results  from  the 
analytical  method  are  nearly  equal  to  those  from  the  heat  balance  integral  method  (HBIM), 
though  it  seems  that  the  HBIM  gives  somewhat  lower  values  of  9/6s  for  smaller  values  of  t) 
and  larger  values  when  r?  is  greater  than  0.5.  A  similar  trend  exists  for  t  =  4.0  (t'  =  1 .0)  ex¬ 
cept  that  the  discrepancy  between  analytical  and  HBIM  grows  larger  and  the  inverse  occurs  at 
much  larger  values  of  ij.  A  similar  statement  can  be  made  for  t  =  20  (t'  =  5.0).  The  HBIM 
gives  much  lower  values  of  0/0„  and  the  inverse  occurs  at  much  higher  values  of 


ALTERNATIVE  METHOD 


As  shown  in  Figure  4  for  t  =  0.4,  4  and  20,  the  deviations  of  HBIM  from  the  analytical 
method  grow  greater  and  greater  as  t  increases.  With  this  type  of  approximation  method, 
one  can  only  expect  a  fair  degree  of  accuracy  for  the  variable  itself,  which  in  this  case  is 


I  ^O)  0( y)  dy. 


Furthermore,  it  is  well  known  that  the  process  of  differentiation  only  increases  the  error. 
Therefore  a  fair  approximation  of  V does  not  necessarily  guarantee  a  similar  fair  approxima¬ 
tion  of  its  derivatives;  in  this  case,  it  means  the  error  for  9  will  be  greater  since 


1  3V 


F(r, i)  dri  ‘ 


Consequently  one  would  expect  a  better  result  if  the  approximation  of  9  can  be  used  instead. 


■  /■ 


Assuming  a  temperature  profile  of 


and  integrating  eq  6  for  0  to  e, 


l  d2d  ,  \  dd  , 

—  =  l  Tr  dr'~ 


F(  n) 


The  final  result  is 


p-  d-e-0 


l_  dt 
3  dr 


(36) 


or 


dt 


6(1 -e  O  ‘ 


(37) 


However,  for  F(rj)  =  e\  an  alternative  solution  with  a  closed  form  can  be  obtained  as 

(37a) 


2  2e* 

t  =  In  f  -  —  +  —  —  Ei(t)  +  >  -  2 


which  is  briefly  derived  in  Appendix  C.  Similarly  for  a  cubic  representation  of  temperature, 
i.e. 


the  dimensionless  temperature  is  found  to  be 


(38) 


The  final  result  is 


(39) 


r 


24(e"'  +  f  -  1) 


(40) 


For  this  case  and  under  the  same  condition,  i.e.  F(rj) 
lution  is 


e’’,  the  alternative  closed  form  so- 


r  =  4(e'-l)+  —  (e'-4)  +  lne-£/(f)  +  7  -  2.5.  (40a) 

f  ( 

Equations  37,  37a,  40  and  40a  are  derived  under  the  condition  that  F(i})  =  e’',  F(0)  =  F'(0) 
=  1 .  Somewhat  detailed  derivations  of  these  two  equations  are  given  in  Appendices  B  and  C. 
Numerical  values  of  r  and  t  evaluated  from  eq  37  and  40  are  listed  in  Table  3  and  graphically 
shown  in  Figure  5  to  facilitate  the  determination  of  the  specific  values  of  r.  For  r  =  0.4,  4 


Table  3.  Values  of  e  at  specific 
values  of  r  from  eq  37  and  40. 


T 

Equation  37 

Equation 

0.01 

0.328 

0.466 

0.04 

— 

— 

0.05 

0.693 

— 

0.08 

— 

1.220 

0.10 

0.943 

1.347 

0.20 

— 

1.822 

0.40 

— 

2.440 

0.50 

1.850 

— 

0.60 

— 

2.880 

0.80 

— 

3.234 

1.0 

2.426 

3.533 

2.0 

3.147 

4.623 

4.0 

4.045 

6.000 

6.0 

4.668 

6.965 

8.0 

5.161 

7.732 

10.0 

5.575 

8.380 

14.0 

6.258 

— 

20.0 

7.067 

10.724 

30.0 

— 

12.366 

40.0 

— 

13.672 

60.0 

— 

15.738 

80.0 

— 

17.381 

1000 

12.148 

18.768 

200  0 

15.316 

23.796 

400.0 

19.303 

30.126 

600.0 

22.099 

34.565 

8  OOO 

24.324 

38.099 

1000.0 

26.203 

41.082 

Figure  5.  r  vs  t  from  numerical  integration 
of  eq  37  and  40. 


and  20,  values  of  e  are  evaluated  from  Figure  5  and  from  eq  37a  and  40a,  and  substituted  in¬ 
to  eq  36  and  39.  The  dimensionless  temperature  ratio  0/0,  is  plotted  in  Figure  4.  It  is  evident 
that,  for  small  values  of  t,  the  results  from  HBIM,  alternative  HBIM  (AHBIM)  and  the  analyt¬ 
ical  are  close  to  each  other.  At  smaller  values  of  rj,  AHBIM  gives  almost  the  same  values  of 
0/0,  as  the  analytical  method.  It  is  also  noted  that  the  simpler  temperature  representation 
(i.e.  eq  36)  gives  a  better  approximation  to  the  analytical  results.  In  general,  AHBIM  always 
gives  higher  ratios  of  0/0,  for  all  q.  As  in  the  HBIM,  the  discrepancy  between  the  AHBIM  and 
the  analytical  results  increases  as  r  increases.  However,  it  can  be  observed  from  Figure  4b 
and  c  that  the  AHBIM  gives  a  much  superior  approximation  than  HBIM  for  smaller  i?  values 
in  all  the  three  cases  considered.  But  it  is  evident  that  for  some  special  functions  of  F{ t?), 
where  a  closed  form  solution  is  feasible,  the  results  are  nearly  identical  to  those  from  the  ana¬ 
lytical  solution  especially  for  large  r  values. 


CONCLUSION  AND  COMMENTS 


The  HBIM  has  been  extended  to  the  solution  of  heat  transfer  problems  in  media  with  varia¬ 
ble  properties,  such  as  snow.  The  essence  of  the  mathematical  treatment  involves  the  intro¬ 
duction  of  a  new  variable  V  and  subsequently  obtaining  the  temperature  profile  by  differen¬ 
tiation.  In  AHBIM,  integration  was  carried  out  directly  and  the  results  thus  obtained  have 


proven  to  be  closer  to  the  analytical  results  than  those  from  the  HB1M.  The  results  obtained 
from  second-degree  polynomial  temperature  representation  were  found  to  be  nearly  as  good 
as  the  results  from  the  cubic  temperature  representation,  though  the  latter  involves  much 
work  in  deriving  eq  4  and  its  numerical  evaluation  of  t.  However,  the  AHBIM  is  more  re¬ 
stricted  because  the  integration  cannot  be  carried  out  unless  the  functional  dependence  of 
F(rj)  is  explicitly  given.  For  a  special  case  of  F(it)  =  e\  a  closed  form  solution  was  obtained, 
which  provides  results  nearly  identical  to  those  from  analytical  methods. 

Only  the  case  of  a  semi-infinite  solid  with  constant  boundary  conditions  has  been  consid¬ 
ered,  but  it  can  be  easily  extended  to  the  case  of  a  finite  slab.  In  such  a  case  the  finite  slab 
would  behave  exactly  the  same  as  a  semi-infinite  solid  until  b  =  L  (L  is  the  thickness  of  the 
slab).  The  conditions  at  v  =  L  would  have  to  be  replaced,  one  of  which  would  be  the  addi¬ 
tional  boundary  condition  at  v  =  1.  and  the  other  determined  from  the  heat  balance  integral. 
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2  +  F  '(0). 


2  +  F'  (0)f 


dr  [  3 [2  +  F 1  (0)t) 


which  can  be  simplified  to 

d<  4  dt 

~  /'( 0)  F  ( 0)  |2  +  F’ (0)<]1 2  • 


Integrating  eq  A9  from  i  =  0  and  <  =  0  to  »,  we  have 
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2  +  F'  (0)t 


IF'  (0)r. 


APPENDIX  B:  DERIV  ATIONS  OF  EQUATIONS  37  AND  40 


Integrating  eq  6  from  0  to  f,  we  have 


\  1  tVd  ,  \  df)  , 

r,  n  ,  dr]  =  —  dr) 

0  r(<l)dr)-  dr) 


if  the  temperature  profile  is  represented  by  eq  36,  i.e. 

. ('  -7 

With  the  use  of  Leibniz’s  rule  for  a  general  function, 


d  */*  .  db  da 

Jr  i  f(xJ)dx  =  fihJ>  Jr  -f(a'l)  Jr 

ait) 

The  right-hand  side  of  eq  6  then  becomes 


j-  \  a  dt)  -  |  y  d>)  *  »(> ) 

drn  .tdr  dr 


<y 

J  /  V' 


lit 


(x.r)  dx. 


but  H(i )  =  0;  therefore 


i  »«.  £r)»a,. 

u  o 


Bv  substitution  of  eq  36  into  B2  and  integration,  eq  B2  becomes 
d 


;w 

:h 


d>) 


The  left-hand  side  of  eq  6  is  integrated  by  parts  with  the  function  /-'(<])  = 


1  ,V"  v  2S,  -s 

.  r  :  ,  -  dr)  =  —  (1  -e  ). 

/  (’))  dr)-  i- 


Combining  eq  B3  and  B4,  i.e. 

2", 


(1  e 


and  Integrating  respectively  for  r  from  0  to  r  and  >  from  0  to  (,  we  have 
V  > 1  d> 

T  =  ’  6(TO  • 


(6) 

(36) 


(Bl) 

(B2) 

(B3) 
( ’’  (see  eq  29)  as 
(B4) 

(B5) 

(37) 


For  the  derivation  of  eq  40,  the  only  difference  is  the  use  of  a  cubic  temperature  profile  in 
stead  of  eq  36,  i.e. 


The  right-hand  side  of  eq  6  becomes 


<M  .  „  dt 

Tn  lhl  ~  ^\4, 


and  the  left-hand  side  is 


'r  1  d‘6  , 

j  F(v)  dv>  dv  ~  t 


Combining  eq  B6  and  B7,  it  follows  that 


Integrating  eq  B8  from  0  to  t  and  0  to  e  respectively, 


=  i  1 


24  *  ((-l+ef) 
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APPENDIX  C:  DERIVATION  OE  EQUATIONS  37a  and  40a 


Integrating  eq  6  from  0  to  <,  i.e. 


1  d:>)  .  ,  , 

I  77,-.  -7-r-  dr)  =  1  -r-  d»j. 

0  t{n)  dry  o  ilrj 


((■>) 


For  a  temperature  profile  of 


*-*.(■-  t)' 


(36) 


and  the  application  of  Leibniz's  rule  (and  from  the  fact  »(>)  =  0)  the  right-hand  side  of  eq  6 
becomes 


]s*-  i\ »*• 


(Cl) 


Equation  6  can  then  be  rewritten  as 


lUd'- 


(C2) 


For  6  given  by  eq  36  and  £(? j)  =  e",  eq  C2  can  be  transformed  to 

?e'\  * 

e  dt  \  f  f !  +  f 2  /  dr  ' 

Upon  differentiation  and  subsequent  integration  eq  C3  becomes 
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ce 


£<(n) 


<C4) 


The  limit  of  the  last  term  can  be  written  as 


lim  [-£/(6)  +  ln6  +  7l  +  Iim  4- (e6- 1) 

A—0  A— 0  o 


(C  5) 


The  first  part  of  eq  C  5  is  zero  and  the  second  part  approaches  a  value  ol  2  when  <i  -  0. 
Therefore,  the  final  expression  for  r  is 


t  =  lnr  -  —  +  —  -£/(<)-*->-  2 

*  f 


(37a) 


where  £/(»)  is  an  exponential  integral  and  is  the  I  uler’s  constant  (-,  0.57722) 

For  the  cubic  temperature  profile 


Following  exactly  the  same  procedures  involved  in  deriving  eq  37a,  eq  6  becomes 
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which  can  be  differentiated  and  then  integrated  to 
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The  last  term  can  be  grouped  as 


lim  [In 0  -  Et(b)  +  7]  +  lim  -  -j-  (e  -4)  . 

0  4_o  0  b 


(C8; 


7  he  first  term  in  eq  C8  is  zero  The  limiting  value  of  the  second  term  can  be  obtained  by  Tay 
lor  expansion  of 
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=  2.5. 


Therefore,  for  the  cubic  temperature  profile,  the  final  expression  of  t  is 
r  =  A  (e  l)  1  (e'-4)  +  lm- £7(0  +  7  -  2.5. 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC  format  is  repro¬ 
duced  below. 
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